home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PIBSIGS / APINVB.PAS < prev    next >
Pascal/Delphi Source File  |  1986-06-22  |  5KB  |  178 lines
  1. C*****APINVB
  2.       DOUBLE PRECISION FUNCTION APINVB( P, Alpha, Beta )
  3. C
  4.       DOUBLE PRECISION       P
  5.       DOUBLE PRECISION       Alpha
  6.       DOUBLE PRECISION       Beta
  7. C
  8. C        Function:  APINVB
  9. C
  10. C        Purpose:  CALCULATES INVERSE Beta DISTRIBUTION
  11. C
  12. C        Calling Sequence:
  13. C
  14. C             B := APINVB( P, Alpha, Beta )
  15. C
  16. C                  P      --- PROBABILITY FOR WHICH PERCENTAGE
  17. C                             POINT IS TO BE FOUND (DOUBLE PRECISION,INP
  18. C                  Alpha  --- FIRST PARAMETER OF Beta DISTRIBUTION
  19. C                             (DOUBLE PRECISION,INPUT)
  20. C                  Beta   --- SECOND PARAMETER OF Beta DISTRIBUTION
  21. C                             (DOUBLE PRECISION,INPUT)
  22. C
  23. C                  B      --- RESULTING PERCENTAGE POINT OF Beta DIST.
  24. C
  25. C        Calls:  CDBeta  (CUMULATIVE Beta DISTRIBUTION)
  26. C
  27. C        Called By:  MANY
  28. C
  29. C        Remarks:
  30. C
  31. C             BECAUSE OF THE RELATIONSHIP BETWEEN THE Beta
  32. C             DISTRIBUTION AND THE F DISTRIBUTION, THIS
  33. C             ROUTINE CAN BE USED TO FIND THE INVERSE OF
  34. C             THE F DISTRIBUTION.
  35. C
  36. C             THE PERCENTAGE POINT IS RETURNED AS -1.0
  37. C             IF ANY ERROR OCCURS.
  38. C
  39. C             THIS ROUTINE USES NEWTONS METHOD
  40. C             TO FIND THE PERCENTAGE POINT.
  41. C
  42. C             THE ALGORITHM IS BASED UPON AS110
  43. C             FROM APPLIED STATISTICS.
  44. C
  45.       INTEGER    ITER,   MAXITR
  46. C
  47.       DOUBLE PRECISION       EPSZ
  48. C
  49.       DOUBLE PRECISION       XIM1,   XI,     XIP1,   FIM1,   FI,     W
  50.       DOUBLE PRECISION       CMPLBT, ADJ,    SQ,     R,      S,      T
  51.       DOUBLE PRECISION       H,      PP,     AA,     BB,     A1,     B1
  52.       DOUBLE PRECISION       G
  53. C
  54.       DOUBLE PRECISION       ALGama
  55.       DOUBLE PRECISION       CDBeta
  56. C
  57.       DOUBLE PRECISION APINVN
  58.       DOUBLE PRECISION PSING
  59. C
  60.       DATA       EPSZ  / 1.0D-26 /
  61. C
  62.       DATA       MAXITR/ 50 /
  63. C
  64. C        (1)  CHECK VALIDITY OF PARAMETERS.
  65. C
  66.       IF( Alpha <= 0.0 ) GOTO 9005
  67.       IF( Beta  <= 0.0 ) GOTO 9005
  68. C
  69.       IF( ( P > 1.0 ) OR ( P < 0.0 ) ) GOTO 9005
  70. C
  71. C        (2)  FLIP PARAMETERS SO THAT
  72. C             'P' FOR EVALUATION IS <:= .5
  73. C
  74.       IF( P > 0.5 ) GOTO 10
  75. C
  76.          AA     := Alpha
  77.          BB     := Beta
  78.          PP     := P
  79.          GOTO 20
  80. C
  81.   10     AA     := Beta
  82.          BB     := Alpha
  83.          PP     := 1.0 - P
  84. C
  85. C        (3)  GENERATE INITIAL APPROXIMATION.
  86. C             SEVERAL DIFFERENT ONES USED, DEPENDING UPON
  87. C             VALUES OF (P,Alpha,Beta).
  88. C
  89.   20  CMPLBT := ALGama( AA ) + ALGama( BB ) -
  90.      1         ALGama( AA + BB )
  91. C
  92.       PSING  := 1.0 - PP
  93.       FI     := APINVN( PSING )
  94. C
  95.       IF( ( AA > 1.0 ) AND ( BB > 1.0 ) ) GOTO 40
  96. C
  97.          R      := BB + BB
  98.          T      := 1.0 / ( 9.0 * BB )
  99.          T      := R * ( 1.0 - T + FI * SQRT( T ) ) ** 3
  100. C
  101.          IF( T > 0.0 ) GOTO 30
  102. C
  103.             XI     := 1.0 - EXP( ( LOG( ( 1.0 - PP ) * BB )
  104.      1                            + CMPLBT ) / BB )
  105.             GOTO 50
  106. C
  107.   30     T      := ( 4.0 * AA + R - 2.0 ) / T
  108. C
  109.          IF( T <= 1.0 ) XI := EXP( (LOG( PP * AA ) + CMPLBT) / PP)
  110.          IF( T > 1.0 ) XI := 1.0 - 2.0 / ( T + 1.0 )
  111.          GOTO 50
  112. C
  113.   40     R      := ( FI * FI - 3.0 ) / 6.0
  114.          S      := 1.0 / ( AA + AA - 1.0 )
  115.          T      := 1.0 / ( BB + BB - 1.0 )
  116.          H      := 2.0 / ( S + T )
  117.          W      := FI * SQRT( H + R ) / H - ( T - S ) *
  118.      1            ( R + 5.0 / 6.0 - 2.0 / ( 3.0 * H ) )
  119.          XI     := AA / ( AA + BB * EXP( W + W ) )
  120. C
  121. C        (4)  BEGIN NEWTON-RAPHSON LOOP.
  122. C
  123.   50  A1     := 1.0 - AA
  124.       B1     := 1.0 - BB
  125.       FIM1   := 0.0
  126.       SQ     := 1.0
  127.       XIM1   := 1.0
  128. C
  129.       DO 100 ITER := 1 , MAXITR
  130. C
  131.          FI     := CDBeta( XI, AA, BB )
  132.          IF( FI < 0.0 ) GOTO 9005
  133. C
  134.          FI     := ( FI - PP ) * EXP( CMPLBT + A1 * LOG( XI ) +
  135.      1                               B1 * LOG( 1.0 - XI ) )
  136. C
  137.          IF( ( FI * FIM1 ) <= 0.0 ) XIM1 := SQ
  138. C
  139.          G      := 1.0
  140. C
  141.   60     ADJ    := G * FI
  142.          SQ     := ADJ * ADJ
  143. C
  144.          IF( SQ < XIM1 ) GOTO 70
  145. C
  146.             G      := G / 3.0
  147.             GOTO 60
  148. C
  149.   70        XIP1     := XI - ADJ
  150.             IF( ( XIP1 .GE. 0.0 ) AND
  151.      1          ( XIP1 <= 1.0 ) ) GOTO 80
  152. C
  153.                G      := G / 3.0
  154.                GOTO 60
  155. C
  156.   80     IF( XIM1    <= EPSZ ) GOTO 9000
  157.          IF( FI * FI <= EPSZ ) GOTO 9000
  158. C
  159.          IF( ( XIP1 = 0.0 ) OR
  160.      1       ( XIP1 = 1.0 ) ) GOTO 70
  161. C
  162.          IF( XIP1 = XI ) GOTO 9000
  163. C
  164.             XI     := XIP1
  165.             FIM1   := FI
  166. C
  167.   100 CONTINUE
  168. C
  169.  9000 APINVB := XI
  170.       IF( P > 0.5 ) APINVB := 1.0 - APINVB
  171. C
  172.       RETURN
  173. C
  174.  9005 XI     := -1.0
  175.       GOTO 9000
  176. C
  177.       END
  178.